#include "DetectorConstruction.hh"
#include "Messenger.hh"
#include "CathodeSD.hh"
#include "PrimaryGeneratorAction.hh"

#include <geant4/G4Material.hh>
#include <geant4/G4Box.hh>
#include <geant4/G4Tubs.hh>
#include <geant4/G4Cons.hh>
#include <geant4/G4Ellipsoid.hh>
#include <geant4/G4SubtractionSolid.hh>
#include <geant4/G4UnionSolid.hh>
#include <geant4/G4LogicalVolume.hh>
#include <geant4/G4PVPlacement.hh>
#include <geant4/G4PVParameterised.hh>
#include <geant4/G4SDManager.hh>
#include <geant4/G4GeometryTolerance.hh>
#include <geant4/G4GeometryManager.hh>
#include <geant4/G4NistManager.hh>
#include <geant4/G4RunManager.hh>
#include <geant4/G4PhysicalVolumeStore.hh>
#include <geant4/G4LogicalVolumeStore.hh>
#include <geant4/G4SolidStore.hh>
#include <geant4/G4LogicalSkinSurface.hh>
#include <geant4/G4LogicalBorderSurface.hh>
#include <geant4/G4SurfaceProperty.hh>
#include <geant4/G4UImanager.hh>
#include <geant4/G4OpticalSurface.hh>
#include <geant4/G4UserLimits.hh>
#include <geant4/G4VisAttributes.hh>
#include <geant4/G4Colour.hh>
#include <geant4/G4ios.hh>
#include <geant4/G4GeneralParticleSource.hh>
#include <geant4/G4RegionStore.hh>

#include <string>
#include <sstream>
#include <limits>

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
DetectorConstruction::DetectorConstruction()
:solidWorld(0),  logicWorld(0),  physiWorld(0),
 solidConicalScint(0), solidEllipsoidalScint(0),
 logicScint(0), physiScint(0),
 solidConicalShellShaft(0), solidConicalShellCap(0),
 solidEllipsoidalShellBase(0), solidEllipsoidalShellInterior(0),solidShell(0),
 logicShell(0),physiShell(0),
 solidAirCap(0),logicAirCap(0),physiAirCap(0),
 solidGrease(0),logicGrease(0),physiGrease(0),
 solidWindow(0), solidShaft(0), solidEnvelope(0),logicEnvelope(0),physiEnvelope(0),
 solidCathode(0),logicCathode(0),physiCathode(0),
 messenger(0), sScintillationMaterial("BaF"),
 sScintillatorFaceSurface("groundAir"), sScintillatorLateralSurface("groundAir"),
 sScintillatorShape("cylindric"),
 fWorldLength(0.), fPMTLength(70.0*mm), fCathodeThickness(1.0*mm),
 fSourceDistance(0.), fScintBaseDiameter(0.), fScintThickness(0.),
 fPMTDiameter(0.), fWindowThickness(0.), fEnvelopeThickness(0.),
 fGreaseThickness(0.), fShellThickness(0.), fAirGap(0.),
  theCathodeSD(0)
{
  messenger = new Messenger(this);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
DetectorConstruction::~DetectorConstruction()
{
  delete messenger;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4VPhysicalVolume* DetectorConstruction::Construct()
{
  G4GeometryManager::GetInstance()->SetWorldMaximumExtent(10*m);
  //G4cout << "Computed tolerance = "
  //       << G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()/mm
  //       << " mm" << G4endl;

  setDefaultDimensions();
  defineMaterials();
  defineSurfaces();
  G4VPhysicalVolume* detector = constructDetector();
  update();
  return detector;
}

void DetectorConstruction::setScintillationMaterial(G4String s)
{
  sScintillationMaterial = s;
}

void DetectorConstruction::setScintillatorFaceSurface(G4String s)
{
  sScintillatorFaceSurface = s;
}

void DetectorConstruction::setScintillatorLateralSurface(G4String s)
{
  sScintillatorLateralSurface = s;
}

void DetectorConstruction::setScintillatorShape(G4String s)
{
  sScintillatorShape = s;
  if (s != "conical")
    fScintHeadDiameter = fScintBaseDiameter;
}

void DetectorConstruction::setSourceDistance(G4double v)
{
  fSourceDistance = v;
}

void DetectorConstruction::setScintDiameter(G4double v)
{
  fScintBaseDiameter = v;
  if (sScintillatorShape != "conical")
    fScintHeadDiameter = v;
}

void DetectorConstruction::setScintHeadDiameter(G4double v)
{
  fScintHeadDiameter = v;
  if (sScintillatorShape != "conical")
    fScintBaseDiameter = v;
}

void DetectorConstruction::setScintThickness(G4double v)
{
  fScintThickness = v;
}

void DetectorConstruction::setPMTDiameter(G4double v)
{
  fPMTDiameter = v;
}

void DetectorConstruction::setWindowThickness(G4double v)
{
  fWindowThickness = v;
}

void DetectorConstruction::setEnvelopeThickness(G4double v)
{
  fEnvelopeThickness = v;
}

void DetectorConstruction::setGreaseThickness(G4double v)
{
  fGreaseThickness = v;
}

void DetectorConstruction::setShellThickness(G4double v)
{
  fShellThickness = v;
}

void DetectorConstruction::setAirGap(G4double v)
{
  fAirGap = v;
}

void DetectorConstruction::setDefaultDimensions()
{
  fSourceDistance = 3.0*cm;
  fScintBaseDiameter = 40*mm;
  fScintHeadDiameter = 40*mm;
  fScintThickness = 12.*mm;
  fPMTDiameter = 50.8*mm;
  fWindowThickness = 2.0*mm;
  fEnvelopeThickness = 2.0*mm;
  fGreaseThickness = 0.1*mm;
  fShellThickness = 0.5*mm;
  fAirGap = 0.2*mm;

  updateWorldLength();
}

void DetectorConstruction::updateDetector()
{
  G4LogicalBorderSurface::CleanSurfaceTable();
  update();
}

void DetectorConstruction::setFileName(G4String fn)
{
  filename = fn;
}

void DetectorConstruction::setSaveOption(CathodeSD::SaveData type, bool save)
{
  if (theCathodeSD)
    theCathodeSD->setSaveOption(type, save);
}

G4String DetectorConstruction::getFileName() const
{
  return filename;
}

void DetectorConstruction::Print(std::ostream &stream)
{
  const PrimaryGeneratorAction *pgenaction = dynamic_cast<const PrimaryGeneratorAction*>(G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
  G4GeneralParticleSource *source = 0;
  if (pgenaction)
    source = pgenaction->getSource();

  stream << "# scintillationMaterial: " << logicScint->GetMaterial()->GetName() << std::endl;
  if (logicScint->GetSolid() != solidEllipsoidalScint)
    stream << "# scintillatorFaceSurface: " << G4LogicalBorderSurface::GetSurface(physiScint, physiAirCap)->GetSurfaceProperty()->GetName() << std::endl;
  stream << "# scintillatorLateralSurface: " << G4LogicalBorderSurface::GetSurface(physiScint, physiWorld)->GetSurfaceProperty()->GetName() << std::endl;
  stream << "# scintillatorShape: ";
  if (logicScint->GetSolid() == solidEllipsoidalScint)
    stream << "ellipsoidal";
  else if (solidConicalScint->GetOuterRadiusMinusZ()==solidConicalScint->GetOuterRadiusPlusZ())
    stream << "cylindric";
  else
    stream << "conical";
  stream << std::endl;
  stream << "# scintDiameter: " << solidConicalScint->GetOuterRadiusPlusZ()*2.0/mm << " mm" << std::endl;
  stream << "# scintHeadDiameter: " << solidConicalScint->GetOuterRadiusMinusZ()*2.0/mm << " mm" << std::endl;
  stream << "# scintThickness: " << solidConicalScint->GetZHalfLength()*2.0/mm << " mm" << std::endl;
  if (source) {
    stream << "# sourceDistance: " << -source->GetCurrentSource()->GetPosDist()->GetCentreCoords().getZ()/mm << " mm" << std::endl;
    stream << "# sourceMinTheta: " << source->GetCurrentSource()->GetAngDist()->GetMinTheta()/degree << "°" << std::endl;
    stream << "# sourceMaxTheta: " << source->GetCurrentSource()->GetAngDist()->GetMaxTheta()/degree << "°" << std::endl;
  }
  stream << "# PMTDiameter: " << solidShaft->GetOuterRadius()*2.0/mm << " mm" << std::endl;
  stream << "# greaseThickness: " << solidGrease->GetZHalfLength()*2.0/mm << " mm" << std::endl;
  stream << "# envelopeThickness: " << (solidShaft->GetOuterRadius() - solidShaft->GetInnerRadius())/mm << " mm" << std::endl;
  stream << "# shellThickness: " << solidConicalShellCap->GetZHalfLength()*2.0/mm << " mm" << std::endl;
  stream << "# windowThickness: " << solidWindow->GetZHalfLength()*2.0/mm << " mm" << std::endl;
  stream << "# airGap: " << solidAirCap->GetZHalfLength()*2.0/mm << " mm" << std::endl;
}

void DetectorConstruction::update()
{
  G4GeometryManager::GetInstance()->OpenGeometry();

  bool isellipsoidal = sScintillatorShape == "ellipsoidal";

  G4ThreeVector positionAirCap = G4ThreeVector(0, 0, fShellThickness+fAirGap/2.0);
  G4ThreeVector positionScint = G4ThreeVector(0, 0, fScintThickness/2.0+fAirGap+fShellThickness);
  G4double shellLength = fScintThickness + fShellThickness + fAirGap;
  G4ThreeVector positionShell = G4ThreeVector(0, 0, fShellThickness/2.0); // boolean solid positioned by cap's center
  if (isellipsoidal) {
    positionScint += G4ThreeVector(0., 0., fScintThickness/2.0);
    positionShell = positionScint;
  }
  G4ThreeVector positionGrease = G4ThreeVector(0,0,fScintThickness+fShellThickness+fAirGap+fGreaseThickness/2.0);
  G4ThreeVector positionCathode = G4ThreeVector(0, 0, fScintThickness+fShellThickness+fAirGap+fGreaseThickness+fWindowThickness+fCathodeThickness/2.0);
  G4ThreeVector positionEnvelope = G4ThreeVector(0, 0, fScintThickness+fShellThickness+fAirGap+fGreaseThickness+fWindowThickness/2.0);

  updateWorldLength();

  // compute values for cylindric/conical shape
  // half apex angle
  G4double alpha = std::atan((fScintBaseDiameter-fScintHeadDiameter)/2.0/fScintThickness);
  G4double airgapxsize = fAirGap / std::cos(alpha);
  G4double shellxsize = fShellThickness / std::cos(alpha);

  G4double scintRadiusBot = fScintBaseDiameter/2.0; //
  G4double scintRadiusTop = fScintHeadDiameter/2.0; //
  G4double airCapRadiusBot = scintRadiusTop; //
  G4double airCapRadiusTop = airCapRadiusBot - fAirGap*std::tan(alpha);
  G4double capRadiusBot = airCapRadiusTop + airgapxsize;
  G4double capRadiusTop = capRadiusBot - fShellThickness*std::tan(alpha);
  G4double shellInnerRadiusBot = scintRadiusBot + airgapxsize; //
  G4double shellInnerRadiusTop = capRadiusTop;
  G4double shellOuterRadiusBot = shellInnerRadiusBot + shellxsize; //
  G4double shellOuterRadiusTop = capRadiusTop + shellxsize;

  assert(alpha == 0.0 || sScintillatorShape == "conical");

  // update dimensions
  solidWorld->SetXHalfLength(0.5*fWorldLength);
  solidWorld->SetYHalfLength(0.5*fWorldLength);
  solidWorld->SetZHalfLength(0.5*fWorldLength);

  solidConicalScint->SetZHalfLength(fScintThickness/2.0);
  solidConicalScint->SetOuterRadiusMinusZ(scintRadiusTop);
  solidConicalScint->SetOuterRadiusPlusZ(scintRadiusBot);
  solidEllipsoidalScint->SetSemiAxis(scintRadiusBot, scintRadiusBot, fScintThickness);
  solidEllipsoidalScint->SetZCuts(-fScintThickness, 0.);
  if (isellipsoidal)
    logicScint->SetSolid(solidEllipsoidalScint);
  else
    logicScint->SetSolid(solidConicalScint);
  physiScint->SetTranslation(positionScint);

  delete solidShell;
  solidConicalShellShaft->SetZHalfLength(shellLength/2.0);
  solidConicalShellShaft->SetInnerRadiusMinusZ(shellInnerRadiusTop);
  solidConicalShellShaft->SetInnerRadiusPlusZ(shellInnerRadiusBot);
  solidConicalShellShaft->SetOuterRadiusMinusZ(shellOuterRadiusTop);
  solidConicalShellShaft->SetOuterRadiusPlusZ(shellOuterRadiusBot);
  solidConicalShellCap->SetZHalfLength(fShellThickness/2.0);
  solidConicalShellCap->SetOuterRadiusMinusZ(capRadiusTop + fShellThickness/2.0); // added half shell tickness to avoid common surface with shaft (overlap is required)
  solidConicalShellCap->SetOuterRadiusPlusZ(capRadiusBot + fShellThickness/2.0);  // added half shell tickness to avoid common surface with shaft (overlap is required)
  solidEllipsoidalShellBase->SetSemiAxis(shellOuterRadiusBot, shellOuterRadiusBot, shellLength);
  solidEllipsoidalShellBase->SetZCuts(-shellLength, 0.);
  solidEllipsoidalShellInterior->SetSemiAxis(shellInnerRadiusBot, shellInnerRadiusBot, shellLength-fShellThickness);
  solidEllipsoidalShellInterior->SetZCuts(-(shellLength-fShellThickness), 0.);
  if (isellipsoidal)
    solidShell = new G4SubtractionSolid("sshell", solidEllipsoidalShellBase, solidEllipsoidalShellInterior);
  else
    solidShell = new G4UnionSolid("sshell",
                                  solidConicalShellCap,
                                  solidConicalShellShaft,
                                  new G4RotationMatrix,
                                  G4ThreeVector(0, 0, (shellLength-fShellThickness)/2.0));
  logicShell->SetSolid(solidShell);
  physiShell->SetTranslation(positionShell);

  solidAirCap->SetZHalfLength(fAirGap/2.0);
  solidAirCap->SetOuterRadiusMinusZ(airCapRadiusTop);
  solidAirCap->SetOuterRadiusPlusZ(airCapRadiusBot);
  delete physiAirCap;
  if (isellipsoidal)
    physiAirCap = 0;
  else
    physiAirCap = new G4PVPlacement(0, positionAirCap, logicAirCap, "pAirCap", logicWorld, false, 0);

  solidGrease->SetZHalfLength(fGreaseThickness/2.0);
  solidGrease->SetOuterRadius(fScintBaseDiameter/2.0);
  physiGrease->SetTranslation(positionGrease);

  delete solidEnvelope;
  solidWindow->SetZHalfLength(fWindowThickness/2.0);
  solidWindow->SetOuterRadius(fPMTDiameter/2.0);
  solidShaft->SetZHalfLength(fPMTLength/2.0);
  solidShaft->SetOuterRadius(fPMTDiameter/2.0);
  solidShaft->SetInnerRadius(fPMTDiameter/2.0 - fEnvelopeThickness);
  solidEnvelope = new G4UnionSolid("senvelope",
                                   solidWindow,
                                   solidShaft,
                                   new G4RotationMatrix,
                                   G4ThreeVector(0, 0, (fPMTLength+fWindowThickness)/2.0)
                                   );
  logicEnvelope->SetSolid(solidEnvelope);
  physiEnvelope->SetTranslation(positionEnvelope);

  solidCathode->SetZHalfLength(fCathodeThickness/2.0);
  solidCathode->SetOuterRadius(fPMTDiameter/2.0 - fEnvelopeThickness);
  physiCathode->SetTranslation(positionCathode);

  //--------- Update gamma source settings
  G4UImanager* UI = G4UImanager::GetUIpointer();
  std::ostringstream sourcepos;
  sourcepos << "/gps/pos/centre 0 0 -" << fSourceDistance/cm << " cm" << std::flush;
  UI->ApplyCommand(sourcepos.str());
  //G4cout << sourcepos << G4endl;
  std::ostringstream maxthetastr;
  G4double maxtheta = std::max(
                               std::atan(shellOuterRadiusTop/fSourceDistance),
                               std::atan(shellOuterRadiusBot/(fSourceDistance+fShellThickness+fAirGap+fScintThickness))
                              );
  maxthetastr << "/gps/ang/maxtheta " << maxtheta << " rad" << std::flush;
  UI->ApplyCommand(maxthetastr.str());
  //G4cout << maxtheta.str() << G4endl;

  // set scintillation material
  logicScint->SetMaterial(G4Material::GetMaterial(sScintillationMaterial));

  // setup surfaces
  if (!isellipsoidal) {
    G4LogicalBorderSurface* scint2AirFaceBorderSurface = 0;
    scint2AirFaceBorderSurface = new G4LogicalBorderSurface("scint2AirFaceBorderSurface", physiScint, physiAirCap, getSurfaceObjectByName(sScintillatorFaceSurface));
  }

  G4LogicalBorderSurface* scint2AirLateralBorderSurface = 0;
  scint2AirLateralBorderSurface = new G4LogicalBorderSurface("scint2AirLateralBorderSurface", physiScint, physiWorld, getSurfaceObjectByName(sScintillatorLateralSurface));

  // shell / cap (not configureable)
  if (!isellipsoidal) {
    G4LogicalBorderSurface* airCap2ShellSurface = 0;
    airCap2ShellSurface = new G4LogicalBorderSurface("air2CapSurface1", physiAirCap, physiShell, airGroundAluminum);
  }
  G4LogicalBorderSurface* air2ShellSurface = 0;
  air2ShellSurface = new G4LogicalBorderSurface("air2ShellSurface", physiWorld, physiShell, airGroundAluminum);

  G4LogicalBorderSurface* envelope2CathodeSurface = 0;
  envelope2CathodeSurface = new G4LogicalBorderSurface("envelope2CathodeSurface", physiEnvelope, physiCathode, silicaCathodeMaterial);

  G4RunManager::GetRunManager()->DefineWorldVolume(physiWorld, true);
  G4RunManager::GetRunManager()->GeometryHasBeenModified();
  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
  G4RegionStore::GetInstance()->UpdateMaterialList(physiWorld);

#ifdef G4VIS_USE
  // update visualisation
  UI->ApplyCommand("/vis/viewer/rebuild");
#endif
}

void DetectorConstruction::defineMaterials()
{
  //--------- Material definition ---------

  G4NistManager* man = G4NistManager::Instance();

  G4Element *C = man->FindOrBuildElement("C");
  G4Element *H = man->FindOrBuildElement("H");
  G4Element *Si = man->FindOrBuildElement("Si");
  G4Element *O = man->FindOrBuildElement("O");
  G4Element *Sb = man->FindOrBuildElement("Sb");
  G4Element *Rb = man->FindOrBuildElement("Rb");
  G4Element *Cs = man->FindOrBuildElement("Cs");
  G4Element *Lu = man->FindOrBuildElement("Lu");
  G4Element *Y = man->FindOrBuildElement("Y");
  G4Element *Ce = man->FindOrBuildElement("Ce");
  G4Element *La = man->FindOrBuildElement("La");
  G4Element *Br = man->FindOrBuildElement("Br");
  G4Element *Na = man->FindOrBuildElement("Na");
  G4Element *I = man->FindOrBuildElement("I");
  G4Element *Tl = man->FindOrBuildElement("Tl");

  //------------------------------
  // Air
  G4Material* Air = man->FindOrBuildMaterial("G4_AIR");
  Air->SetName("Air");
  const G4int numentries = 2;
  G4double energies[numentries] = { 0.1*eV, 10.0*eV };
  G4double vacrindices[numentries] = { 1., 1. };
  G4double airabsorpti[numentries] = { 100.*m, 100.*m }; // avoid infinite light-paths
  G4MaterialPropertiesTable* airprop = new G4MaterialPropertiesTable();
  airprop->AddProperty("ABSLENGTH", energies, airabsorpti, numentries);
  airprop->AddProperty("RINDEX", energies, vacrindices, numentries);
  Air->SetMaterialPropertiesTable(airprop);
  //------------------------------

  //------------------------------
  // Aluminum
  G4Material* Aluminum = man->FindOrBuildMaterial("G4_Al");
  Aluminum->SetName("Aluminum");
  //------------------------------

  //------------------------------
  // common BaF2
  const G4int numentriesbaf = 10;
  G4double bafenergies[numentriesbaf] = { 1.2*eV, 3.99*eV, 4.00*eV, 4.01*eV, 4.02*eV, 5.63*eV, 5.64*eV, 6.36*eV, 6.37*eV, 6.5*eV }; // saint-gobain
  G4double baffastcomp[numentriesbaf] = { 0.0,    0.0,     0.0,     0.0,     0.0,     0.0,     0.67,    0.33,    0.0,     0.0    }; // saint-gobain
  G4double bafslowcomp[numentriesbaf] = { 0.0,    0.0,     1.0,     1.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.0    }; // saint-gobain
  G4double bafrindices[numentriesbaf] = { 1.50,   1.50,    1.50,    1.50,    1.50,    1.54,    1.54,    1.54,    1.54,    1.54   }; // saint-gobain
  G4double bafabsorpti[numentriesbaf] = { 34*cm,  34*cm,   34*cm,   34*cm,   34*cm,   34*cm,   34*cm,   34*cm,   34*cm,   34*cm  }; // Ma1992-2
  G4Material *BaFtemplate = man->FindOrBuildMaterial("G4_BARIUM_FLUORIDE");

  // Perfect BaF2
  G4Material* PerfBaF = new G4Material("PerfectBaF", BaFtemplate->GetDensity(), BaFtemplate, kStateSolid);
  G4MaterialPropertiesTable* perfbafprop = new G4MaterialPropertiesTable;
  perfbafprop->AddProperty("FASTCOMPONENT", bafenergies, baffastcomp, numentriesbaf);
  perfbafprop->AddProperty("SLOWCOMPONENT", bafenergies, bafslowcomp, numentriesbaf);
  perfbafprop->AddProperty("RINDEX",        bafenergies, bafrindices, numentriesbaf);
  perfbafprop->AddProperty("ABSLENGTH",     bafenergies, bafabsorpti,  numentriesbaf);
  perfbafprop->AddConstProperty("SCINTILLATIONYIELD",100/keV); // perfect
  perfbafprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  perfbafprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // perfect
  perfbafprop->AddConstProperty("SLOWTIMECONSTANT",0.000001*ns); // perfect
  perfbafprop->AddConstProperty("YIELDRATIO",0.15); // saint-gobain
  PerfBaF->SetMaterialPropertiesTable(perfbafprop);

  // Perfect Fast BaF2
  G4Material* PerfFastBaF = new G4Material("PerfectFastBaF", BaFtemplate->GetDensity(), BaFtemplate, kStateSolid);
  G4MaterialPropertiesTable* perffastbafprop = new G4MaterialPropertiesTable;
  perffastbafprop->AddProperty("FASTCOMPONENT", bafenergies, baffastcomp, numentriesbaf);
  perffastbafprop->AddProperty("RINDEX",        bafenergies, bafrindices, numentriesbaf);
  perffastbafprop->AddProperty("ABSLENGTH",     bafenergies, bafabsorpti,  numentriesbaf);
  perffastbafprop->AddConstProperty("SCINTILLATIONYIELD",100/keV); // perfect
  perffastbafprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  perffastbafprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // perfect
  perffastbafprop->AddConstProperty("YIELDRATIO",1.0); // only fast fraction
  PerfFastBaF->SetMaterialPropertiesTable(perffastbafprop);

  // Immediate BaF2
  G4Material* ImmediateBaF = new G4Material("ImmediateBaF", BaFtemplate->GetDensity(), BaFtemplate, kStateSolid);
  G4MaterialPropertiesTable* immbafprop = new G4MaterialPropertiesTable;
  immbafprop->AddProperty("FASTCOMPONENT", bafenergies, baffastcomp, numentriesbaf);
  immbafprop->AddProperty("SLOWCOMPONENT", bafenergies, bafslowcomp, numentriesbaf);
  immbafprop->AddProperty("RINDEX",        bafenergies, bafrindices, numentriesbaf);
  immbafprop->AddProperty("ABSLENGTH",     bafenergies, bafabsorpti,  numentriesbaf);
  immbafprop->AddConstProperty("SCINTILLATIONYIELD",11.8/keV); // saint-gobain
  immbafprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  immbafprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // immediate
  immbafprop->AddConstProperty("SLOWTIMECONSTANT",0.000001*ns); // immediate
  immbafprop->AddConstProperty("YIELDRATIO",0.15); // saint-gobain
  ImmediateBaF->SetMaterialPropertiesTable(immbafprop);

  // Immediate Fast BaF2
  G4Material* ImmediateFastBaF = new G4Material("ImmediateFastBaF", BaFtemplate->GetDensity(), BaFtemplate, kStateSolid);
  G4MaterialPropertiesTable* immfastbafprop = new G4MaterialPropertiesTable;
  immfastbafprop->AddProperty("FASTCOMPONENT", bafenergies, baffastcomp, numentriesbaf);
  immfastbafprop->AddProperty("RINDEX",        bafenergies, bafrindices, numentriesbaf);
  immfastbafprop->AddProperty("ABSLENGTH",     bafenergies, bafabsorpti,  numentriesbaf);
  immfastbafprop->AddConstProperty("SCINTILLATIONYIELD",1.77/keV); // saint-gobain * fastyield
  immfastbafprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  immfastbafprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // immediate
  immfastbafprop->AddConstProperty("YIELDRATIO",1.0); // only fast fraction
  ImmediateFastBaF->SetMaterialPropertiesTable(immfastbafprop);

  // BaF2
  G4Material* BaF = new G4Material("BaF", BaFtemplate->GetDensity(), BaFtemplate, kStateSolid);
  G4MaterialPropertiesTable* bafprop = new G4MaterialPropertiesTable;
  bafprop->AddProperty("FASTCOMPONENT", bafenergies, baffastcomp, numentriesbaf);
  bafprop->AddProperty("SLOWCOMPONENT", bafenergies, bafslowcomp, numentriesbaf);
  bafprop->AddProperty("RINDEX",        bafenergies, bafrindices, numentriesbaf);
  bafprop->AddProperty("ABSLENGTH",     bafenergies, bafabsorpti,  numentriesbaf);
  bafprop->AddConstProperty("SCINTILLATIONYIELD",11.8/keV); // saint-gobain
  bafprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  bafprop->AddConstProperty("FASTTIMECONSTANT",0.7*ns); // saint-gobain
  bafprop->AddConstProperty("SLOWTIMECONSTANT",630.*ns); // saint-gobain
  bafprop->AddConstProperty("YIELDRATIO",0.15); // saint-gobain
  BaF->SetMaterialPropertiesTable(bafprop);

  // Fast BaF2
  G4Material* FastBaF = new G4Material("FastBaF", BaFtemplate->GetDensity(), BaFtemplate, kStateSolid);
  G4MaterialPropertiesTable* fastbafprop = new G4MaterialPropertiesTable;
  fastbafprop->AddProperty("FASTCOMPONENT", bafenergies, baffastcomp, numentriesbaf);
  fastbafprop->AddProperty("RINDEX",        bafenergies, bafrindices, numentriesbaf);
  fastbafprop->AddProperty("ABSLENGTH",     bafenergies, bafabsorpti,  numentriesbaf);
  fastbafprop->AddConstProperty("SCINTILLATIONYIELD",1.77/keV); // saint-gobain * fastyield
  fastbafprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  fastbafprop->AddConstProperty("FASTTIMECONSTANT",0.7*ns); // saint-gobain
  fastbafprop->AddConstProperty("YIELDRATIO",1.0); // only fast fraction
  FastBaF->SetMaterialPropertiesTable(fastbafprop);
  //------------------------------

  //------------------------------
  // common LYSO
  const G4int numentrieslyso = 6;
  G4double lysoenergies[numentrieslyso] = { 1.2*eV, 2.94*eV, 2.95*eV, 2.96*eV, 2.97*eV, 6.5*eV }; // saint-gobain (420 nm)
  G4double lysofastcomp[numentrieslyso] = { 0.0, 0., 1.0, 1.0, 0., 0.0 };
  G4double lysorindices[numentrieslyso] = { 1.81, 1.81, 1.81, 1.81, 1.81, 1.81 }; // saint-gobain
  G4double lysoabsorptionlength[numentrieslyso] = { 42*cm, 42*cm, 42*cm, 42*cm, 42*cm, 42*cm }; // vilardi2006
  G4Material* LYSOtemplate = new G4Material("LYSOtemplate", 7.1*g/cm3, 5, kStateSolid);
  LYSOtemplate->AddElement(Lu, 71.43*perCent);
  LYSOtemplate->AddElement(Y, 4.03*perCent);
  LYSOtemplate->AddElement(Si, 6.37*perCent);
  LYSOtemplate->AddElement(O, 18.14*perCent);
  LYSOtemplate->AddElement(Ce, 0.02*perCent); // cooke2000

  // Perfect LYSO
  G4Material* PerfLYSO = new G4Material("PerfectLYSO", LYSOtemplate->GetDensity(), LYSOtemplate, kStateSolid);
  G4MaterialPropertiesTable* perflysoprop = new G4MaterialPropertiesTable();
  perflysoprop->AddProperty("FASTCOMPONENT", lysoenergies, lysofastcomp, numentrieslyso);
  perflysoprop->AddProperty("RINDEX",        lysoenergies, lysorindices, numentrieslyso);
  perflysoprop->AddProperty("ABSLENGTH",     lysoenergies, lysoabsorptionlength,  numentrieslyso);
  perflysoprop->AddConstProperty("SCINTILLATIONYIELD",100./keV); // perfect
  perflysoprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  perflysoprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // perfect
  perflysoprop->AddConstProperty("YIELDRATIO",1.0);
  PerfLYSO->SetMaterialPropertiesTable(perflysoprop);

  // Immediate LYSO
  G4Material* ImmediateLYSO = new G4Material("ImmediateLYSO", LYSOtemplate->GetDensity(), LYSOtemplate, kStateSolid);
  G4MaterialPropertiesTable* immlysoprop = new G4MaterialPropertiesTable();
  immlysoprop->AddProperty("FASTCOMPONENT", lysoenergies, lysofastcomp, numentrieslyso);
  immlysoprop->AddProperty("RINDEX",        lysoenergies, lysorindices, numentrieslyso);
  immlysoprop->AddProperty("ABSLENGTH",     lysoenergies, lysoabsorptionlength,  numentrieslyso);
  immlysoprop->AddConstProperty("SCINTILLATIONYIELD",32./keV); // saint-gobain
  immlysoprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  immlysoprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // immediate
  immlysoprop->AddConstProperty("YIELDRATIO",1.0);
  ImmediateLYSO->SetMaterialPropertiesTable(immlysoprop);

  // LYSO
  G4Material* LYSO = new G4Material("LYSO", LYSOtemplate->GetDensity(), LYSOtemplate, kStateSolid);
  G4MaterialPropertiesTable* lysoprop = new G4MaterialPropertiesTable();
  lysoprop->AddProperty("FASTCOMPONENT", lysoenergies, lysofastcomp, numentrieslyso);
  lysoprop->AddProperty("RINDEX",        lysoenergies, lysorindices, numentrieslyso);
  lysoprop->AddProperty("ABSLENGTH",     lysoenergies, lysoabsorptionlength,  numentrieslyso);
  lysoprop->AddConstProperty("SCINTILLATIONYIELD",32./keV); // saint-gobain
  lysoprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  lysoprop->AddConstProperty("FASTTIMECONSTANT",41.0*ns); // saint-gobain
  lysoprop->AddConstProperty("YIELDRATIO",1.0);
  LYSO->SetMaterialPropertiesTable(lysoprop);
  //------------------------------

  //------------------------------
  // common LaBr
  const G4int numentrieslabr = 6;
  G4double labrenergies[numentrieslabr] = { 1.2*eV, 3.25*eV, 3.26*eV, 3.27*eV, 3.28*eV, 6.5*eV }; // saint-gobain (420 nm)
  G4double labrfastcomp[numentrieslabr] = { 0.0, 0., 1.0, 1.0, 0., 0.0 };
  G4double labrrindices[numentrieslabr] = { 1.9, 1.9, 1.9, 1.9, 1.9, 1.9 }; // saint-gobain
  G4double labrabsorptionlength[numentrieslabr] = { 1.0*m, 1.0*m, 1.0*m, 1.0*m, 1.0*m, 1.0*m }; // guessed
  G4Material* LaBrtemplate = new G4Material("LaBrtemplate", 5.08*g/cm3, 3, kStateSolid); // saint-gobain
  LaBrtemplate->AddElement(La, 36.62*perCent);
  LaBrtemplate->AddElement(Br, 63.19*perCent);
  LaBrtemplate->AddElement(Ce, 0.19*perCent); // loef2001

  // Perfect LaBr
  G4Material* PerfLaBr = new G4Material("PerfectLaBr", LaBrtemplate->GetDensity(), LaBrtemplate, kStateSolid);
  G4MaterialPropertiesTable* perflabrprop = new G4MaterialPropertiesTable();
  perflabrprop->AddProperty("FASTCOMPONENT", labrenergies, labrfastcomp, numentrieslabr);
  perflabrprop->AddProperty("RINDEX",        labrenergies, labrrindices, numentrieslabr);
  perflabrprop->AddProperty("ABSLENGTH",     labrenergies, labrabsorptionlength,  numentrieslabr);
  perflabrprop->AddConstProperty("SCINTILLATIONYIELD",100./keV); // perfect
  perflabrprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  perflabrprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // perfect
  perflabrprop->AddConstProperty("YIELDRATIO",1.0);
  PerfLaBr->SetMaterialPropertiesTable(perflabrprop);

  // Immediate LaBr
  G4Material* ImmediateLaBr = new G4Material("ImmediateLaBr", LaBrtemplate->GetDensity(), LaBrtemplate, kStateSolid);
  G4MaterialPropertiesTable* immlabrprop = new G4MaterialPropertiesTable();
  immlabrprop->AddProperty("FASTCOMPONENT", labrenergies, labrfastcomp, numentrieslabr);
  immlabrprop->AddProperty("RINDEX",        labrenergies, labrrindices, numentrieslabr);
  immlabrprop->AddProperty("ABSLENGTH",     labrenergies, labrabsorptionlength,  numentrieslabr);
  immlabrprop->AddConstProperty("SCINTILLATIONYIELD",75.0/keV); // dehaas2008
  immlabrprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  immlabrprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // immediate
  immlabrprop->AddConstProperty("YIELDRATIO",1.0);
  ImmediateLaBr->SetMaterialPropertiesTable(immlabrprop);

  // LaBr
  G4Material* LaBr = new G4Material("LaBr", LaBrtemplate->GetDensity(), LaBrtemplate, kStateSolid);
  G4MaterialPropertiesTable* labrprop = new G4MaterialPropertiesTable();
  labrprop->AddProperty("FASTCOMPONENT", labrenergies, labrfastcomp, numentrieslabr);
  labrprop->AddProperty("RINDEX",        labrenergies, labrrindices, numentrieslabr);
  labrprop->AddProperty("ABSLENGTH",     labrenergies, labrabsorptionlength,  numentrieslabr);
  labrprop->AddConstProperty("SCINTILLATIONYIELD",75.0/keV); // dehaas2008
  labrprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  labrprop->AddConstProperty("FASTTIMECONSTANT",16.0*ns); // saint-gobain
  labrprop->AddConstProperty("YIELDRATIO",1.0);
  LaBr->SetMaterialPropertiesTable(labrprop);
  //------------------------------

  //------------------------------
  // common NaI
  const G4int numentriesnai = 6;
  G4double naienergies[numentriesnai] = { 1.2*eV, 2.98*eV, 2.99*eV, 3.0*eV, 3.01*eV, 6.5*eV }; // saint-gobain (415 nm)
  G4double naifastcomp[numentriesnai] = { 0.0, 0.0, 1.0, 1.0, 0.0, 0.0 };
  G4double naislowcomp[numentriesnai] = { 0.0, 0.0, 1.0, 1.0, 0.0, 0.0 };
  G4double nairindices[numentriesnai] = { 1.85, 1.85, 1.85, 1.85, 1.85, 1.85 }; // saint-gobain
  G4double naiabsorptionlength[numentriesnai] = { 1.0*m, 1.0*m, 1.0*m, 1.0*m, 1.0*m, 1.0*m }; // guessed
  G4Material* NaItemplate = new G4Material("NaItemplate", 3.67*g/cm3, 3, kStateSolid); // saint-gobain
  NaItemplate->AddElement(Na, 15.32*perCent);
  NaItemplate->AddElement(I, 84.59*perCent);
  NaItemplate->AddElement(Tl, 0.08*perCent); // trefilova2002

  // Perfect NaI
  G4Material* PerfNaI = new G4Material("PerfectNaI", NaItemplate->GetDensity(), NaItemplate, kStateSolid);
  G4MaterialPropertiesTable* perfnaiprop = new G4MaterialPropertiesTable();
  perfnaiprop->AddProperty("FASTCOMPONENT", naienergies, naifastcomp, numentriesnai);
  perfnaiprop->AddProperty("RINDEX",        naienergies, nairindices, numentriesnai);
  perfnaiprop->AddProperty("ABSLENGTH",     naienergies, naiabsorptionlength,  numentriesnai);
  perfnaiprop->AddConstProperty("SCINTILLATIONYIELD",100./keV); // perfect
  perfnaiprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  perfnaiprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // perfect
  perfnaiprop->AddConstProperty("YIELDRATIO",1.0);
  PerfNaI->SetMaterialPropertiesTable(perfnaiprop);

  // Immediate NaI
  G4Material* ImmediateNaI = new G4Material("ImmediateNaI", NaItemplate->GetDensity(), NaItemplate, kStateSolid);
  G4MaterialPropertiesTable* immnaiprop = new G4MaterialPropertiesTable();
  immnaiprop->AddProperty("FASTCOMPONENT", naienergies, naifastcomp, numentriesnai);
  immnaiprop->AddProperty("RINDEX",        naienergies, nairindices, numentriesnai);
  immnaiprop->AddProperty("ABSLENGTH",     naienergies, naiabsorptionlength,  numentriesnai);
  immnaiprop->AddConstProperty("SCINTILLATIONYIELD",38.0/keV); // saint-gobain
  immnaiprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  immnaiprop->AddConstProperty("FASTTIMECONSTANT",0.000001*ns); // immediate
  immnaiprop->AddConstProperty("YIELDRATIO",1.0);
  ImmediateNaI->SetMaterialPropertiesTable(immnaiprop);

  // NaI
  G4Material* NaI = new G4Material("NaI", NaItemplate->GetDensity(), NaItemplate, kStateSolid);
  G4MaterialPropertiesTable* naiprop = new G4MaterialPropertiesTable();
  naiprop->AddProperty("FASTCOMPONENT", naienergies, naifastcomp, numentriesnai);
  naiprop->AddProperty("SLOWCOMPONENT", naienergies, naislowcomp, numentriesnai);
  naiprop->AddProperty("RINDEX",        naienergies, nairindices, numentriesnai);
  naiprop->AddProperty("ABSLENGTH",     naienergies, naiabsorptionlength,  numentriesnai);
  naiprop->AddConstProperty("SCINTILLATIONYIELD",38.0/keV); // saint-gobain
  naiprop->AddConstProperty("RESOLUTIONSCALE",1.0);
  naiprop->AddConstProperty("FASTTIMECONSTANT",250.0*ns); // saint-gobain
  naiprop->AddConstProperty("SLOWTIMECONSTANT",250.0*ns); // saint-gobain
  naiprop->AddConstProperty("SLOWSCINTILLATIONRISETIME",27.0*ns); // weber2000
  naiprop->AddConstProperty("YIELDRATIO",0.5);
  NaI->SetMaterialPropertiesTable(naiprop);
  //------------------------------

  //------------------------------------------------
  // Polydimethylsiloxane (Grease)
  G4Material* Polydimethylsiloxane = new G4Material("Polydimethylsiloxane", 0.97*g/cm3, 4, kStateLiquid);
  Polydimethylsiloxane->AddElement(Si, 1);
  Polydimethylsiloxane->AddElement(O, 1);
  Polydimethylsiloxane->AddElement(C, 2);
  Polydimethylsiloxane->AddElement(H, 6);
  G4MaterialPropertiesTable* polydimethylsiloxaneprop = new G4MaterialPropertiesTable();
  const G4int numentriespolydimethylsiloxane = 3;
  G4double polydimethylsiloxaneenergy[numentriespolydimethylsiloxane] = {1.2*eV, 3.1*eV, 6.5*eV};
  G4double polydimethylsiloxaneabsorp[numentriespolydimethylsiloxane] = {10.*cm, 10.*cm, 10.*cm};
  G4double polydimethylsiloxanerindex[numentriespolydimethylsiloxane] = {1.4, 1.4, 1.4};
  polydimethylsiloxaneprop->AddProperty("ABSLENGTH", polydimethylsiloxaneenergy, polydimethylsiloxaneabsorp, numentriespolydimethylsiloxane);
  polydimethylsiloxaneprop->AddProperty("RINDEX", polydimethylsiloxaneenergy, polydimethylsiloxanerindex, numentriespolydimethylsiloxane);
  Polydimethylsiloxane->SetMaterialPropertiesTable(polydimethylsiloxaneprop);
  //------------------------------------------------

  //------------------------------
  // Fused silica
  G4Material* FusedSilica = new G4Material("FusedSilica", 2.201*g/cm3, 2, kStateSolid);
  FusedSilica->AddElement(Si, 1);
  FusedSilica->AddElement(O, 2);
  G4MaterialPropertiesTable* fusedsilicaprop = new G4MaterialPropertiesTable();
  const G4int numentriesfusedsilica = 3;
  G4double fusedsilicaenergy[numentriesfusedsilica] = {1.2*eV, 3.1*eV, 6.5*eV};
  G4double fusedsilicaabsorp[numentriesfusedsilica] = {2.*m, 2.*m, 2.*m};
  G4double fusedsilicarindex[numentriesfusedsilica] = {1.56, 1.47, 1.45};
  fusedsilicaprop->AddProperty("ABSLENGTH", fusedsilicaenergy, fusedsilicaabsorp, numentriesfusedsilica);
  fusedsilicaprop->AddProperty("RINDEX", fusedsilicaenergy, fusedsilicarindex, numentriesfusedsilica);
  FusedSilica->SetMaterialPropertiesTable(fusedsilicaprop);
  //------------------------------

  //------------------------------
  // Bialkali Cathode (dummy)
  G4Material* BialkaliCathode = new G4Material("BialkaliCathode", 3*g/cm3, 3, kStateSolid);
  BialkaliCathode->AddElement(Sb, 1);
  BialkaliCathode->AddElement(Rb, 1);
  BialkaliCathode->AddElement(Cs, 1);
  G4MaterialPropertiesTable* bialkalicathodeprop = new G4MaterialPropertiesTable();
//  const G4int numentriesbialkalicath = 2;
//  G4double bialkalicathodeenergy[numentriesbialkalicath] = {1.2*eV, 6.5*eV};
//  G4double bialkalicathodeabsorp[numentriesbialkalicath] = {1.e-6*mm, 1.e-6*mm}; // absorb all
//  bialkalicathodeprop->AddProperty("ABSLENGTH", bialkalicathodeenergy, bialkalicathodeabsorp, numentriesbialkalicath);
  bialkalicathodeprop->AddProperty("RINDEX", fusedsilicaenergy, fusedsilicarindex, numentriesfusedsilica); // use values from window to prevent refraction
  BialkaliCathode->SetMaterialPropertiesTable(bialkalicathodeprop);
  //------------------------------

  //G4cout << *(G4Material::GetMaterialTable()) << G4endl;
}


void DetectorConstruction::defineSurfaces()
{
  G4double ener[2] = {.1*eV, 10.*eV};

  BGOPolishedAirTeflon = new G4OpticalSurface("BGOPolishedAirTeflon", unified);
  BGOPolishedAirTeflon->SetType(dielectric_LUT);
  BGOPolishedAirTeflon->SetModel(LUT);
  BGOPolishedAirTeflon->SetFinish(polishedteflonair);

  BGOGroundAirTeflon = new G4OpticalSurface("BGOGroundAirTeflon", unified);
  BGOGroundAirTeflon->SetType(dielectric_LUT);
  BGOGroundAirTeflon->SetModel(LUT);
  BGOGroundAirTeflon->SetFinish(groundteflonair);

  polishedAir = new G4OpticalSurface("polishedAir", unified);
  polishedAir->SetType(dielectric_dielectric);
  polishedAir->SetModel(unified);
  polishedAir->SetFinish(ground); // necessary even for polished surfaces to enable UNIFIED code
  polishedAir->SetSigmaAlpha(1.3 * degree); // Janecek2010


  groundAir = new G4OpticalSurface("groundAir", unified);
  groundAir->SetType(dielectric_dielectric);
  groundAir->SetModel(unified);
  groundAir->SetFinish(ground);
  groundAir->SetSigmaAlpha(0.2 * radian); // Janecek2010

  // common paint attributes
  G4double paintspecularlobe[2] = {1.0, 1.0};
  G4double paintspecularspike[2] = {0.0, 0.0};
  G4double paintbackscatter[2] = {0.0, 0.0};
  G4double paintrindex[2] = {1.61, 1.61}; // Janecek2010

  // white painted
  G4double whitepaintrefl[2] = {0.955, 0.955}; // Janecek2010

  polishedWhitePainted = new G4OpticalSurface("polishedWhitePainted", unified);
  polishedWhitePainted->SetType(dielectric_dielectric);
  polishedWhitePainted->SetModel(unified);
  polishedWhitePainted->SetFinish(groundbackpainted);
  polishedWhitePainted->SetSigmaAlpha(1.3 * degree); // Janecek2010
  G4MaterialPropertiesTable* polishedWhitePaintedProperty = new G4MaterialPropertiesTable();
  polishedWhitePaintedProperty->AddProperty("RINDEX",ener,paintrindex,2);
  polishedWhitePaintedProperty->AddProperty("SPECULARLOBECONSTANT",ener,paintspecularlobe,2);
  polishedWhitePaintedProperty->AddProperty("SPECULARSPIKECONSTANT",ener,paintspecularspike,2);
  polishedWhitePaintedProperty->AddProperty("BACKSCATTERCONSTANT",ener,paintbackscatter,2);
  polishedWhitePaintedProperty->AddProperty("REFLECTIVITY",ener,whitepaintrefl,2);
  polishedWhitePainted->SetMaterialPropertiesTable(polishedWhitePaintedProperty);

  groundWhitePainted = new G4OpticalSurface("groundWhitePainted", unified);
  groundWhitePainted->SetType(dielectric_dielectric);
  groundWhitePainted->SetModel(unified);
  groundWhitePainted->SetFinish(groundbackpainted);
  groundWhitePainted->SetSigmaAlpha(12.0 * degree); // Janecek2010
  G4MaterialPropertiesTable* groundWhitePaintedProperty = new G4MaterialPropertiesTable();
  groundWhitePaintedProperty->AddProperty("RINDEX",ener,paintrindex,2);
  groundWhitePaintedProperty->AddProperty("SPECULARLOBECONSTANT",ener,paintspecularlobe,2);
  groundWhitePaintedProperty->AddProperty("SPECULARSPIKECONSTANT",ener,paintspecularspike,2);
  groundWhitePaintedProperty->AddProperty("BACKSCATTERCONSTANT",ener,paintbackscatter,2);
  groundWhitePaintedProperty->AddProperty("REFLECTIVITY",ener,whitepaintrefl,2);
  groundWhitePainted->SetMaterialPropertiesTable(groundWhitePaintedProperty);

  // black painted
  G4double blackpaintrefl[2] = {0.043, 0.043}; // Dury2006

  polishedBlackPainted = new G4OpticalSurface("polishedBlackPainted", unified);
  polishedBlackPainted->SetType(dielectric_dielectric);
  polishedBlackPainted->SetModel(unified);
  polishedBlackPainted->SetFinish(groundbackpainted);
  polishedBlackPainted->SetSigmaAlpha(1.3 * degree); // Janecek2010
  G4MaterialPropertiesTable* polishedBlackPaintedProperty = new G4MaterialPropertiesTable();
  polishedBlackPaintedProperty->AddProperty("RINDEX",ener,paintrindex,2);
  polishedBlackPaintedProperty->AddProperty("SPECULARLOBECONSTANT",ener,paintspecularlobe,2);
  polishedBlackPaintedProperty->AddProperty("SPECULARSPIKECONSTANT",ener,paintspecularspike,2);
  polishedBlackPaintedProperty->AddProperty("BACKSCATTERCONSTANT",ener,paintbackscatter,2);
  polishedBlackPaintedProperty->AddProperty("REFLECTIVITY",ener,blackpaintrefl,2);
  polishedBlackPainted->SetMaterialPropertiesTable(polishedBlackPaintedProperty);

  groundBlackPainted = new G4OpticalSurface("groundBlackPainted", unified);
  groundBlackPainted->SetType(dielectric_dielectric);
  groundBlackPainted->SetModel(unified);
  groundBlackPainted->SetFinish(groundbackpainted);
  groundBlackPainted->SetSigmaAlpha(12.0 * degree); // Janecek2010
  G4MaterialPropertiesTable* groundBlackPaintedProperty = new G4MaterialPropertiesTable();
  groundBlackPaintedProperty->AddProperty("RINDEX",ener,paintrindex,2);
  groundBlackPaintedProperty->AddProperty("SPECULARLOBECONSTANT",ener,paintspecularlobe,2);
  groundBlackPaintedProperty->AddProperty("SPECULARSPIKECONSTANT",ener,paintspecularspike,2);
  groundBlackPaintedProperty->AddProperty("BACKSCATTERCONSTANT",ener,paintbackscatter,2);
  groundBlackPaintedProperty->AddProperty("REFLECTIVITY",ener,blackpaintrefl,2);
  groundBlackPainted->SetMaterialPropertiesTable(groundBlackPaintedProperty);

  // surface definition of the aluminum shell around the scintillation crystal
  airGroundAluminum = new G4OpticalSurface("airGroundAluminum", unified);
  airGroundAluminum->SetType(dielectric_metal);
  airGroundAluminum->SetModel(unified);
  airGroundAluminum->SetFinish(ground);
  airGroundAluminum->SetSigmaAlpha(0.2);
  G4MaterialPropertiesTable *airGroundAluminumProperty = new G4MaterialPropertiesTable();
  G4double refl[2] = {0.7, 0.7};
  airGroundAluminumProperty->AddProperty("REFLECTIVITY",ener,refl,2);
  airGroundAluminum->SetMaterialPropertiesTable(airGroundAluminumProperty);

  silicaCathodeMaterial = new G4OpticalSurface("silicaCathodeMaterial", unified);
  silicaCathodeMaterial->SetType(dielectric_metal);
  silicaCathodeMaterial->SetModel(unified);
  silicaCathodeMaterial->SetFinish(ground);
  silicaCathodeMaterial->SetSigmaAlpha(0.);
  G4MaterialPropertiesTable *silicaCathodeMaterialProperty = new G4MaterialPropertiesTable();
  G4double cathoderefl[2] = {0., 0.};
  G4double cathodeeff[2] = {1., 1.};
  silicaCathodeMaterialProperty->AddProperty("REFLECTIVITY",ener,cathoderefl,2);
  silicaCathodeMaterialProperty->AddProperty("EFFICIENCY",ener,cathodeeff,2);
  silicaCathodeMaterial->SetMaterialPropertiesTable(silicaCathodeMaterialProperty);

}


G4VPhysicalVolume * DetectorConstruction::constructDetector()
{
  G4UImanager* UI = G4UImanager::GetUIpointer();

  UI->ApplyCommand("/gps/pos/type Point");
  UI->ApplyCommand("/gps/ang/type iso");
  UI->ApplyCommand("/gps/ene/mono 66 keV");
  UI->ApplyCommand("/gps/particle gamma");
  UI->ApplyCommand("/gps/ang/rot1 -1 0 0");

  //--------- Determine positions ---------

  G4ThreeVector positionScint = G4ThreeVector(0, 0, fScintThickness/2.0+fShellThickness);
  G4ThreeVector positionGrease = G4ThreeVector(0,0,fScintThickness+fShellThickness+fGreaseThickness/2.0);
  G4ThreeVector positionShell = G4ThreeVector(0, 0, fShellThickness/2.0);
  G4ThreeVector positionCathode = G4ThreeVector(0, 0, fScintThickness+fShellThickness+fGreaseThickness+fWindowThickness+fCathodeThickness/2.0);
  G4ThreeVector positionEnvelope = G4ThreeVector(0, 0, fScintThickness+fShellThickness+fGreaseThickness+fWindowThickness/2.0);
  if (fSourceDistance > fScintBaseDiameter)
    fWorldLength= fSourceDistance*2.0;
  else
    fWorldLength= fScintBaseDiameter*2.0;

  //--------- Definitions of Solids, Logical Volumes, Physical Volumes ---------

  //------------------------------
  // World
  //------------------------------

  G4double HalfWorldLength = 0.5*fWorldLength;
  solidWorld= new G4Box("sworld",HalfWorldLength,HalfWorldLength,HalfWorldLength);
  logicWorld= new G4LogicalVolume( solidWorld, G4Material::GetMaterial("Air"), "lWorld", 0, 0, 0);

  //  Must place the World Physical volume unrotated at (0,0,0).
  //
  physiWorld = new G4PVPlacement(0,               // no rotation
                                 G4ThreeVector(), // at (0,0,0)
                                 logicWorld,      // its logical volume
                                 "pWorld",         // its name
                                 0,               // its mother  volume
                                 false,           // no boolean operations
                                 0);              // copy number

  // Set user cuts to avoid deadlocks
  G4double maxStep = 10.0*m, maxLength = 10.0*m, maxTime = 100.0*ns, minEkin = 0.5*eV;
  logicWorld->SetUserLimits(new G4UserLimits(maxStep,maxLength,maxTime,minEkin));

  //------------------------------
  // Scintillator
  //------------------------------

  solidConicalScint = new G4Cons("sconicalscintillator", 0.*cm, 1.*cm, 0.*cm, 1.*cm, 1.*cm, 0.*deg, 360.*deg);
  solidEllipsoidalScint = new G4Ellipsoid("sellipsoidalscintillator", 1.*cm, 1.*cm, 1.*cm, 0., 0.);
  logicScint = new G4LogicalVolume(solidConicalScint,G4Material::GetMaterial("BaF"),"lScintillator",0,0,0);
  physiScint = new G4PVPlacement(0,               // no rotation
                                  positionScint,  // at (x,y,z)
                                  logicScint,     // its logical volume
                                  "pScintillator",        // its name
                                  logicWorld,      // its mother  volume
                                  false,           // no boolean operations
                                  0);              // copy number

  //------------------------------
  // Shell / Cap
  //------------------------------

  solidConicalShellShaft = new G4Cons("sconicalshell", 0.*cm, 1.*cm, 0.*cm, 1.*cm, 1.*cm, 0.*deg, 360.*deg);
  solidConicalShellCap = new G4Cons("scap", 0.*cm, 1.*cm, 0.*cm, 1.*cm, 1.*cm, 0.*deg, 360.*deg);
  solidEllipsoidalShellBase = new G4Ellipsoid("sellipsoidalshellbase", 2.*cm, 2.*cm, 2.*cm, 0., 0.);
  solidEllipsoidalShellInterior = new G4Ellipsoid("sellipsoidalshellinterior", 1.*cm, 1.*cm, 1.*cm, 0., 0.);
  solidShell = new G4UnionSolid("sshell",
                                solidConicalShellCap,
                                solidConicalShellShaft,
                                new G4RotationMatrix,
                                G4ThreeVector(0, 0, 1*cm));
  logicShell = new G4LogicalVolume(solidShell, G4Material::GetMaterial("Aluminum"), "lShell", 0,0,0);
  physiShell = new G4PVPlacement(0, positionShell, logicShell, "pShell", logicWorld, false, 0);
  solidAirCap = new G4Cons("saircap", 0.*cm, 1.*cm, 0.*cm, 1.*cm, 1.*cm, 0.*deg, 360.*deg);
  logicAirCap = new G4LogicalVolume(solidAirCap, G4Material::GetMaterial("Air"), "lAirCap", 0,0,0);
  //physiAirCap = new G4PVPlacement(0, positionCap, logicAirCap, "pAirCap", logicWorld, false, 0);

  //------------------------------------------------
  // Grease
  //------------------------------------------------

  solidGrease = new G4Tubs("sgrease", 0.*cm, 1.*cm, 1.*cm, 0.*deg, 360.*deg);
  logicGrease = new G4LogicalVolume(solidGrease, G4Material::GetMaterial("Polydimethylsiloxane"), "lGrease", 0,0,0);
  physiGrease = new G4PVPlacement(0,               // no rotation
                                  positionGrease,  // at (x,y,z)
                                  logicGrease,     // its logical volume
                                  "pGrease",        // its name
                                  logicWorld,      // its mother  volume
                                  false,           // no boolean operations
                                  0);              // copy number

  //------------------------------
  // PMT Window / Envelope
  //------------------------------

  solidWindow = new G4Tubs("swindow", 0.*cm, 1.*cm, 1.*cm, 0.*deg, 360.*deg);
  solidShaft = new G4Tubs("sshaft", 0.*cm, 1.*cm, 1.*cm, 0.*deg, 360.*deg);
  solidEnvelope = new G4UnionSolid("senvelope",
                                   solidWindow,
                                   solidShaft,
                                   new G4RotationMatrix,
                                   G4ThreeVector(0, 0, (fPMTLength+fWindowThickness)/2.0)
                                   );
  logicEnvelope = new G4LogicalVolume(solidEnvelope, G4Material::GetMaterial("FusedSilica"), "lEnvelope", 0,0,0);
  physiEnvelope = new G4PVPlacement(0,               // no rotation
                                  positionEnvelope,  // at (x,y,z)
                                  logicEnvelope,     // its logical volume
                                  "pEnvelope",        // its name
                                  logicWorld,      // its mother  volume
                                  false,           // no boolean operations
                                  0);              // copy number

  //------------------------------
  // PMT Cathode
  //------------------------------

  solidCathode = new G4Tubs("scath", 0.*cm, 1.*cm, 1.*cm, 0.*deg, 360.*deg);
  logicCathode = new G4LogicalVolume(solidCathode, G4Material::GetMaterial("BialkaliCathode"), "lCath", 0,0,0);
  physiCathode = new G4PVPlacement(0,               // no rotation
                                  positionCathode,  // at (x,y,z)
                                  logicCathode,     // its logical volume
                                  "pCath",        // its name
                                  logicWorld,      // its mother  volume
                                  false,           // no boolean operations
                                  0);              // copy number

  //------------------------------------------------
  // Define sensitive detector
  //------------------------------------------------

  G4SDManager* SDman = G4SDManager::GetSDMpointer();
  G4String sensitiveDetectorName = "/detector/sensitiveDetector";
  theCathodeSD = new CathodeSD(sensitiveDetectorName, physiCathode);
  SDman->AddNewDetector( theCathodeSD );
  logicEnvelope->SetSensitiveDetector( theCathodeSD );

  //--------- Visualization attributes -------------------------------

  G4VisAttributes* ScintVisAtt = new G4VisAttributes(G4Colour(1.0,0.0,1.0,0.6));
  G4VisAttributes* ShellVisAtt = new G4VisAttributes(G4Colour(0.8,0.8,0.8,0.2));
  G4VisAttributes* CathodeVisAtt = new G4VisAttributes(G4Colour(0.5,0.5,0.5,0.0));
  CathodeVisAtt->SetForceWireframe(true);
  G4VisAttributes* PMTVisAtt = new G4VisAttributes(G4Colour(0.38823529 ,0.67843137, 0.81568627, 0.16));
  G4VisAttributes* GreaseVisAtt = new G4VisAttributes(G4Colour(1.0,1.0,0.0,0.8));
  logicWorld->SetVisAttributes(G4VisAttributes::Invisible);
  logicScint->SetVisAttributes(ScintVisAtt);
  logicShell->SetVisAttributes(ShellVisAtt);
  logicCathode->SetVisAttributes(CathodeVisAtt);
  logicEnvelope->SetVisAttributes(PMTVisAtt);
  logicAirCap->SetVisAttributes(G4VisAttributes::Invisible);
  logicGrease->SetVisAttributes(GreaseVisAtt);

  return physiWorld;
}

G4OpticalSurface * DetectorConstruction::getSurfaceObjectByName(G4String name) const
{
  if (name == "BGOPolishedAirTeflon")
    return BGOPolishedAirTeflon;
  if (name == "BGOGroundAirTeflon")
    return BGOGroundAirTeflon;
  if (name == "polishedAir")
    return polishedAir;
  if (name == "groundAir")
    return groundAir;
  if (name == "polishedWhitePainted")
    return polishedWhitePainted;
  if (name == "groundWhitePainted")
    return groundWhitePainted;
  if (name == "polishedBlackPainted")
    return polishedBlackPainted;
  if (name == "groundBlackPainted")
    return groundBlackPainted;
  return 0;
}

void DetectorConstruction::updateWorldLength()
{
  G4double detectorlength = fScintThickness+fShellThickness+fAirGap+fGreaseThickness+fWindowThickness+fPMTLength;
  fWorldLength = std::max(detectorlength, fSourceDistance)*2.0;
  fWorldLength = std::max(fWorldLength, fPMTDiameter);
  fWorldLength = std::max(fWorldLength, (fScintBaseDiameter/2.0 + fAirGap + fShellThickness)*2.0);
  fWorldLength *= 2.0;
}



